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Sequential  Conjugate  Gradient-Restoration  Algorithm 
for  Optimal  Control  Problems 
with  Nondifferential  Constraints,  Part  1,  Theory^ 

by 

J.R.  CLOUTIER^,  B.P.  MOHANTY^ , and  A.  MIELe'^ 

Abstract.  A sequential  conjugate  gradient-restoration  algorithm 
is  developed  in  order  to  solve  optimal  control  problems  involving 
a functional  subject  to  differential  constraints,  nondifferential 
constraints,  and  terminal  constraints.  The  algorithm  is  composed 
of  a sequence  of  cycles,  each  cycle  consisting  of  two  phases,  a 
conjugate  gradient  phase  and  a restoration  phase. 
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The  conjugate  gradient  phase  involves  a single  iteration 
and  is  designed  to  decrease  the  value  of  the  functional  while 
satisfying  the  constraints  to  first  order.  During  this  iter- 
ation, the  first  variation  of  the  functional  is  minimized,  sub- 
ject to  the  linearized  constraints.  The  minimization  is 
performed  over  the  class  of  variations  of  the  control  and  the 
parameter  which  are  equidistant  from  some  constant  multiple  of 
the  corresponding  variations  of  the  previous  conjugate  gradient 
phase.  For  the  special  case  of  a quadratic  functional  subject 
to  linear  constraints,  various  orthogonality  and  conjugacy 
conditions  hold. 

The  restoration  phase  involves  one  or  more  iterations  and 
is  designed  to  restore  the  constraints  to  a predetermined 
accuracy,  while  the  norm  of  the  variations  of  the  control  and 
the  parameter  is  minimized,  subject  to  the  linearized  constraints. 
The  restoration  phase  is  terminated  whenever  the  norm  of  the 
constraint  error  is  less  than  some  predetermined  tolerance. 

The  sequential  conjugate  gradient-restoration  algorithm  is 
characterized  by  two  main  properties.  First,  at  the  end  of  each 
conjugate  gradient-restoration  cycle,  the  trajectory  satisfies 
the  constraints  to  a given  accuracy;  thus,  a sequence  of  feasible 
suboptimal  solutions  is  produced.  Second,  the  conjugate  gradient 
stepsize  and  the  restoration  stepsize  can  be  chosen  so  that  the 
restoration  phase  preserves  the  descent  property  of  the  conjugate 
gradient  phase;  thus,  the  value  of  the  functional  at  the  end  of 
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any  cycle  is  smaller  than  the  value  of  the  functional  at  the 
beginning  of  that  cycle.  Of  course,  restarting  the  algorithm 
might  be  occasionally  necessary. 

To  facilitate  numerical  integrations,  the  interval  of 
integration  is  normalized  to  unit  length.  Variable-time 
terminal  conditions  are  transformed  into  fixed-time  terminal 
conditions.  Then,  the  actual  time  at  which  the  terminal  bound- 
ary is  reached  becomes  a component  of  a vector  parameter  being 
optimized . 

Convergence  is  attained  whenever  both  the  norm  of  the 
constraint  error  and  the  norm  of  the  error  in  the  optimality 
conditions  are  less  than  some  predetermined  tolerances.  Seve- 
ral numerical  examples  illustrating  the  theory  of  this  paper 
are  given  in  Part  2. 

Key  Words.  Optimal  control,  gradient  methods,  conjugate-gradient 
methods,  numerical  methods,  computing  methods,  gradient-restoration 
algorithms,  sequential  gradient-restoration  algorithms,  sequential 
conjugate  gradient-restoration  algorithms,  nondifferential  cons- 
traints. 
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1 . Introduction 

Approximately  ten  years  ago,  conjugate  gradient  techniques 
began  appearing  on  the  optimal  control  scene.  In  1967,  Lasdon 
et  al  (Ref.  1)  extended  the  conjugate  gradient  method  developed 
by  Fletcher  and  Reeves  for  mathematical  programming  problems  to 
optimal  control  problems.  About  the  same  time,  Horwitz  and 
Sarachik  (Ref.  2)  extended  Davidon's  method  to  a real  Hilbert 
space  and  applied  the  extension  to  a control  problem  with  qua- 
dratic cost  and  linear  constraints.  Shortly  thereafter,  Lasdon 
(Ref.  3)  and  Tripathi  and  Narendra  (Ref.  4)  also  derived  exten- 
sions of  Davidon's  method. 

One  common  limitation  of  the  above  algorithms  is  that  they 
are  not  applicable  directly  to  constrained  control  problems 
(that  is,  problems  involving  terminal  constraints  and/or  bounds 
on  the  state  or  the  control).  However,  these  algorithms  can 
handle  indirectly  constrained  control  problems,  after  conversion 
of  these  problems  to  unconstrained  form;  this  conversion  is 
usually  achieved  by  means  of  penalty  functions. 

Conjugate  gradient  algorithms  which  can  solve  directly 
certain  types  of  constrained  control  problems  were  presented  in 
Refs.  5-8.  Sinnott  and  Luenberger  constructed  an  algorithm  for 
solving  problems  with  linear  terminal  constraints  (Ref.  5) ; 
Heideman  and  Levy  developed  an  algorithm  for  problems  with  ar- 
bitrary terminal  constraints  (Refs.  6-7);  and  Pagurek  and 
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Woodside  constructed  an  algorithm  for  problems  with  bounded 
controls  (Ref.  8) . 

In  the  area  of  ordinary  gradient  methods,  Miele  et  al 
(Ref.  9)  developed  a sequential  gradient-restoration 
algorithm  for  optimal  control  problems  where  the  state  x(t), 
the  control  u(t),  and  the  parameter  tt  must  satisfy  not 
only  differential  constraints  and  terminal  constraints, 
but  also  nondifferential  constraints  everywhere  along  the  in- 
terval of  integration.  The  importance  of  Ref.  9 lies  in  that 
(i)  many  optimization  problems  arise  directly  in  the  form  con- 
sidered there,  (ii)  problems  involving  equality  constraints 
can  be  reduced  to  that  form  through  suitable  transformations, 
and  (iii)  problems  involving  inequality  constraints  can  be  re- 
duced to  that  form  through  suitable  transformations.  Thus,  an 
extremely  large  class  of  problems  can  be  handled.  This  includes 
problems  with  bounded  control,  bounded  state,  bounded  time  rate 
of  change  of  the  state,  as  well  as  problems  where  a bound  is 

imposed  on  some  function  of  the  parameter,  the  control,  the  state, 
and  the  time  rate  of  change  of  the  state. 

This  report  combines  the  ideas  of  Ref.  6 and  those  of  Ref.  9 . 
The  result  is  a sequential  conjugate  gradient-restoration  algo- 
rithm which  can  handle  constrained  minimization  problems,  char- 
acterized by  the  presence  of  nondifferential  constraints,  without 
resorting  to  penalty  functions.  The  algorithm  is  composed  of  a 
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sequence  of  cycles,  each  cycle  consisting  of  two  phases,  a 
conjugate  gradient  phase  and  a restoration  phase. 

The  conjugate  gradient  phase  involves  a single  iteration 
and  is  designed  to  decrease  the  value  of  the  functional  while 
satisfying  the  constraints  to  first  order.  During  this  iter- 
ation, the  first  variation  of  the  functional  is  minimized 
subject  to  the  linearized  constraints.  The  minimization  is 
performed  over  the  class  of  variations  of  the  control  and  the 
parameter  which  are  equidistant  from  some  constant  multiple  of 
the  corresponding  variations  of  the  previous  conjugate  gradient 
phase . 

The  restoration  phase  involves  one  or  more  iterations  and 
is  designed  to  restore  the  constraints  to  a predetermined  accu- 
racy while  the  norm  of  the  variations  of  the  control  and  the 
parameter  is  minimized,  subject  to  the  linearized  constraints. 

The  sequential  conjugate  gradient-restoration  algorithm  is 
characterized  by  two  main  properties.  First,  at  the  end  of 
each  conjugate  gradient-restoration  cycle,  the  trajectory  satis- 
fies the  constraints  to  a given  accuracy;  thus,  a sequence  of 
feasible  suboptimal  solutions  is  produced.  Second,  the  conju- 
gate gradient  stepsize  and  the  restoration  stepsize  can  be 
chosen  so  that  the  restoration  phase  preserves  the  descent 
property  of  the  conjugate  gradient  phase;  thus,  the  value  of 
the  functional  at  the  end  of  any  cycle  is  smaller  than  the  value 
of  the  functional  at  the  beginning  of  that  cycle.  Of  course, 
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restarting  the  algorithm  might  be  occasionally  necessary.  For 
a discussion  of  the  basic  properties  of  the  sequential  gradient- 
restoration  algorithm,  see  Ref.  10. 
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2 . Formulation  of  the  Problem 

Consider  the  problem  of  minimizing  the  functional 


/: 


1=1  f*  (x,u,TT^,0)d0  + [g*  (x,TTj^,0 


(1) 


with  respect  to  the  state  x(G),  the  control  u(0),  and  the  para- 
meters TT*  and  t which  satisfy  the  differential  constraint 


dx/de  = ({)*  (x,u,TT^  ,0),  0<G<t, 


(2) 


the  nondifferential  constraint 


S* (x, u, n* , 9)  = 0 , O<0<T, 


(3) 


and  the  boundary  conditions 


(x)Q  = given,  [ (x , , 0)]^  = 0 . 


(4) 


In  the  above  equations,  the  functions  f*  and  g*  are 
scalar,  the  function  4>^  is  an  n-vector,  the  function  S*  is  a 
k-vector,  and  the  function  is  a q-vector.  The  independent 
variable  is  the  actual  time  0,  while  the  dependent  variables 
are  the  state  x (an  n-vector) , the  control  u (an  m-vector) , 
the  parameter  (a  p*-vector)  , and  the  parameter  t (a  scalar). 
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At  the  initial  time'  0 = 0,  the  n components  of  the  vector  x are 
specified.  At  the  final  time  6=1,  q scalar  relations  are 
specified,  where  q f n + if  t is  fixed  and  q < n + p*  + 1 if  x is  free. 

To  facilitate  the  implementation  of  the  algorithm  on  a 
digital  computer,  we  replace  the  actual  time  6 with  the  nor- 
malized time  t.  The  latter  is  defined  in  such  a way  that  the 
interval  of  integration  has  unit  length.  Thus,  in  normalized 
form,  t=  0 denotes  the  time  at  which  the  initial  boundary  (4-1) 
is  left  and  t = 1 denotes  the  time  at  which  the  terminal  bounda- 
ry (4-2)  is  reached.  The  following  linear  relation  allows  the 
passage  from  the  normalized  time  t to  the  actual  time  0: 

e = xt,  0 < t < 1.  (5) 

The  fact  that  the  normalized  final  time  is  fixed  (t=l) 
does  not  cause  any  loss  of  generality  in  the  problem.  If  the 
actual  final  time  is  free,  x simply  becomes  a parameter  to  be 
optimized  in  the  transformed  problem.  In  view  of  this,  we 
define  the  augmented  parameter  tt  ( a p-vector)  , 

TT  = TT^  or  77  = , (6) 

where  (6-1) holds  if  i is  fixed  and  (6-2)  holds  if  i is  free. 

In  addition  to  the  normalized  time  t and  the  parameter  ti  , 
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we  define  the  following  functions; 


U,7rj^,Tt)  , 

0 < 

t < 

1, 

(7) 

u,  71*  , 1 t)  , 

0 < 

t < 

1, 

(8) 

‘,7T*,Tt), 

0 < 

t < 

1, 

(9) 

(10) 

* , T t)  . 

(11) 

Under  the  above  transformations  and  definitions,  problem 
(l)-(4)  can  be  reformulated  as  follows.  Minimize  the  functional 


'i: 


f {x,u,-n,t)dt  + [g(x,-n,t)] 


(12) 


with  respect  to  the  state  x(t),  the  control  u(t),  and  the  para- 
meter TT  which  satisfy  the  differential  constraint 


X = (}!  {x,U,TT,t)  , 


0 < t < 1 , 


(13) 


the  nondifferential  constraint 


S(X,U,1T,t)  = 0, 


0 < t < 1 , 


(14: 


and  the  boundary  conditions 


[ip(x,v,t)  ]^=  0, 


(x) Q = given, 


(15) 
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From  calculus  of  variations,  v\'e  know  that  the  problem 
(12) -(15)  is  one  of  the  Bolza  type;  it  can  be  recast  as  that 
of  minimizing  the  augmented  functional^ 


= 1 (>  " X + H)  dt  + (G)^ 

•^0 

= (a'^x)^  + { (h-a'^x)  dt  + (G) 

0 Jr, 


with  respect  to  the  state  x(t),  the  control  u(t),  and  the 
parameter  which  satisfy  (13)-(15),  where  the  functions  H and 
G are  given  by 


rp  m 

H = + p^s, 


G = g + n il-, 


and  where  A(t)  is  a variable  Lagrange  multiplier  (an  n-vector) , 
o(t)  is  a variable  Lagrange  multiplier  (a  k-vector)  , and  is 
a constant  Lagrange  multiplier  (a  q-vector) . 

The  optimal  solution  must  satisfy  (13)  — (15)  and  the  first* 
order  optimality  conditions,  namely,  the  Euler  equations 


A = II  , 
X 


II  = 0, 
u 


J 0 


,,dt  + (G^)^  = 0 (181 


In  Eq.  (16),  it  is  tacitly  assumed  that  the  initial  condition 
(15-1)  is  satisfied.  The  second  form  of  Eq . (16)  arises  after 

the  customary  integration  by  parts  is  performed. 
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and  the  following  natural  condition  arising  from  the  transver- 
sality  condition: 


(A  + = 0.  (19) 

Summarizing,  we  seek  functions  x(t),  u(t),  tt  and  multipliers 
A(t),  p(t),  u which  satisfy  the  constraints  (13) -(15)  and  the 
optimality  conditions  (20) -(23): 


A-f  +4>  A-S  p=0, 

X X X 

^u  ■ ‘^u^  SuP= 

lo  +S^p)dt  + (g 

(A+g^  + = 0. 


0 < t < 1 , 

(20) 

0 < t < 1 , 

(21) 

+ 1 = 0, 

(22) 

(23) 

2.1.  Approximate  Methods.  Since  the  differential  system 
(13) -(15)  and  (20) -(23)  is  generally  nonlinear,  some  iterative 
technique  must  be  employed  in  its  solution.  For  this  purpose, 
let  us  define  the  scalar  functionals  P and  Q, which  denote  the 
constraint  error  and  the  error  in  the  optimality  conditions, 
respectively.  We  have 


N (x  - 4) ) dt  + 


1 ' 


N(S)dt  + N(i)^). 


(24) 
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” Jo  £ N(f^  - 0^A  + S^p)dt 

£ 


^^TT  ■ ^ s^p)dt+  (g^+>i^^u)i 


+N  (A+g^+4;^y; 


where  N(b)  denotes  the  norm  squared  of  a vector,  i.e., 


N (b)  = b b 


for  a given  vector  b. 

Note  that,  for  the  optimal  solution,  P=0  and  Q=0. 
an  approximation  to  the  optimal  solution, 

P<Ci,  Q<e2^ 


(25) 


(26) 


For 


(27) 


where  and  are  small,  preselected  numbers. 
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3 . Construction  of  the  Sequential  Conjugate  Gradient- 

Restoration  Algorithm. 

The  sequential  conjugate  gradient-restoration  algorithm 
is  an  iterative  technique  which  includes  a sequence  of  cycles 
having  the  following  properties.^ 

Property  3.1.  The  functions  x{t),  u{t),  tt  available 
both  at  the  beginning  and  at  the  end  of  each  cycle  must  be 
feasible;  that  is,  they  must  be  consistent  with  the  constraints 
(13)-(15)  within  the  preselected  accuracy  (27-1). 

Property  3.2.  The  functions  x{t),  u(t),  ir  produced  at  the 
end  of  each  cycle  must  be  characterized  by  a value  of  the 
functional  I [see  Eg.  (12)]  which  is  smaller  than  that  associ- 
ated with  the  functions  available  at  the  beginning  of  the 
cycle. 

Property  3.3.  The  functions  x(t),  u(t),  tt  produced  at  the 
end  of  each  cycle  must  be  characterized  by  a value  of  the 
augmented  functional  J [see  Eq.  (16)]  which  is  smaller  than 
that  associated  with  the  functions  available  at  the  beginning 
of  the  cycle. 


^Note  that  Property  3.3  is  a consequence  of  Properties  3.1  and 
3.2.  Conversely,  Property  3.2  is  a consequence  of  Properties 
3.1  and  3.3. 
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To  achieve  the  above  properties,  each  cycle  is  made  of 
two  phases,  a conjugate  gradient  phase  and  a restoration  phase. 

Conjugate  Gradient  Phase.  This  phase  is  started  only 
when  the  constraint  error  P satisfies  Ineq.  (27-1).  It  involves 
a single  iteration,  which  is  designed  to  decrease  the  value  of 
the  functional  I or  the  augmented  functional  J,  while  satisfy- 
ing the  constraints  to  first  order.  During  this  iteration,  the 
first  variation  of  the  functional  I is  minimized,  subject  to 
the  linearized  constraints.  The  minimization  is  performed  over 
the  class  of  variations  of  the  control  and  the  parameter  which 
are  equidistant  from  some  constant  multiple  of  the  corresponding 
variations  of  the  previous  conjugate  gradient  phase. 

Restoration  Phase.  This  phase  is  started  only  when  the 
constraint  error  P violates  Ineq.  (27-1).  The  restoration  phase 
involves  one  or  more  iterations.  In  each  restorative  iteration, 
the  objective  is  to  reduce  the  constraint  error  P,  while  the  cons- 
traints are  satisfied  to  first  order  and  the  norm  of  the 
variations  of  the  control  and  the  parameter  is  minimized.  The 
restoration  phase  is  terminated  whenever  Ineq.  (27-1)  is 
satisfied . 

Remark.  During  each  conjugate  gradient  iteration  or  res- 
torative iteration,  use  of  nonlinear  equations  must  be  avoided. 
Therefore,  the  exact  feasibility  equations  (13) -(15)  are  re- 
placed by  their  corresponding  linearized  approximations. These 
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linearized  approximations  do  not  include  forcing  terms  in  the 
conjugate  gradient  phase,  while  they  do  include  forcing  terms 
in  the  restoration  phase. 

Notation . For  any  iteration  of  the  conjugate  gradient 
phase  or  the  restoration  phase,  the  following  terminology  is 
adopted:  x(t),  u(t),  ir  denote  the  nominal  functions;  x(t),  u(t), 

IT  denote  the  varied  functions;  and  Ax(t),  Au(t),  Att  denote  the 
displacements  leading  from  the  nominal  functions  to  the  varied 
functions.  These  quantities  satisfy  the  definitions 


X (t)  = X (t)  + Ax  (t),  u (t)  = u(t)  + Au(t),  TT=7r+A7T.  (28) 

i If  the  variations  appearing  in  (28)  are  linear  in  the  stepsize 

I 

I a,  where  a > 0,  they  take  the  form 

i 

I 

■ Ax  (t)  = aA{t),  Au  (t)  = aB  (t),  Att  = aC,  (29) 

1 

I 

with  the  implication  that 

i 

x(t)=x(t)+aA(t),  u(t)=u(t)+aB(t),  Tr  = TT+aC.  (30) 

The  functions  Ax(t),  Au(t),  Air  must  be  determined  so  as  to 
produce  some  desirable  effect  at  every  iteration,  namely,  the 
decrease  of  the  functionals  I and/or  J and/or  P.  Thus,  the 


L 
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following  descent  properties  are  required: 

i < I , and/or  J < J , and/or  P < P,  (31) 

where  I,J,P  are  associated  with  the  nominal  functions  and  i,J, 

P are  associated  with  the  varied  functions.  Inequalities 
(31-1)  and  (31-2)  characterize  the  conjugate  gradient  phase, 
and  Ineq.  (31-3)  characterizes  the  restoration  phase. 

In  turn,  relations  (31)  can  be  enforced  at  every  itera- 
tion providing  the  stepsize  a is  sufficiently  small  and  the 
functions  A(t),  3(t),  C are  chosen  so  that 

6l  < 0,  and/or  6J  < 0,  and/or  6P  < 0,  (32) 

where  the  symbol  6(...)  denotes  the  first  variation.  Inequali- 
ties (32-1)  and  (32-2)  characterize  the  conjugate  gradient 
phase,  and  Ineq.  (32-3)  characterizes  the  restoration  phase. 

Clearly,  every  iteration  includes  two  basic  operations: 

(a)  the  determination  of  functions  A(t),  B(t),  C consistent  with 
the  first  variation  requirements  (32);  and  (b)  the  determination 
of  the  stepsize  a consistent  with  the  total  variation  require- 
ments ( 31 ) . 

Outline.  In  Section  4,  we  describe  the  equations  of  the 
restoration  phase;  we  show  how  nominal  functions  consistent  with 
the  feasibility  equations  (13) -(15)  can  be  obtained.  In 
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Section  5,  we  describe  the  general  equations  of  the  conjugate 
gradient  phase;  the  linear  case  (case  where  the  constraints 
are  linear)  is  treated  in  Section  6;  the  linear-quadratic  case 
(case  v;here  the  functional  is  quadratic  and  the  constraints 
are  linear)  is  treated  in  Section  7;  here,  we  show  that  cer- 
tain general  conjugacy  and  orthogonality  conditions  hold. 
Always  with  reference  to  the  conjugate  gradient  phase,  the 
nonlinear-nonquadratic  case  (case  where  the  functional  is 
nonquadratic  and/or  the  constraints  are  nonlinear)  is  treated 
in  Section  8;  here,  we  discuss  the  implementation  of  a 
first-order  algorithm  (this  is  an  algorithm  v\?hich  uses  first 
derivatives  at  most),  in  Section  9,  we  discuss  the  descent 
property  of  a complete  conjugate  gradient-restoration  cycle. 

In  Section  10,  we  present  a summary  of  the  sequential  conju- 
gate gradient-restoration  algorithm.  Finally, in  Section  11, 
we  list  the  safeguards  necessary  to  its  implementation  on  a 
digital  computer. 
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4 . Restoration  Phase 

As  stated  before,  the  restoration  phase  is  started  only 
when  the  constraint  error  P violates  Ineq.  (27-1).  The  resto- 
ration phase  involves  one  or  more  iterations.  In  each  restora- 
tive iteration,  the  objective  is  to  reduce  the  constraint  error 
P,  while  the  constraints  are  satisfied  to  first  order  and  the 
norm  of  the  variations  of  the  control  and  the  parameter  is 
minimized.  The  restoration  phase  is  terminated  whenever  Ineq. 
(27-1)  is  satisfied. 

There  are  two  situations  where  the  restoration  phase  is 
employed:  (a)  at  the  very  beginning  of  the  algorithm,  one  needs 
to  generate  nominal  functions  con^'istent  with  the  feasibility 
equations  (13) -(15);  and  (b)  subsequently,  one  needs  to  correct 
possible  constraint  violations  occurring  during  a conjugate 
gradient  phase:  these  constraint  violations  are  due  to  the  fact 
that  Eqs.  (13)- (15)  are  considered  only  in  linearized  form 
during  a conjugate  gradient  phase. 

Linearized  Equations.  Let  x(t),  u(t),  r denote  nominal  func- 
tions not  satisfying  (13)  - (15),  and  let  x(t),  u(t),  tt  denote  varied 
functions  satisf ying ( 1 3 ) - ( 1 5).  To  f irst  order  , the  perturba- 

tions Ax(t),  Au(t),  Att  must  satisfy  the  linearized  constraint 
equations 

Ax  - Ax  - ({)'^Au  - a'^Att  + a {x  - ^)  = 0 , 0<t<l  (33) 

X U IT  _ _ » 
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S Ax  + S Au  + S Att  + cxS  = 0,  0 < t < 1,  (34) 

(Ax)q=0,  (i|;^Ax  + i);^Ait  + onj;)  ^ = 0,  (35) 

where  a denotes  a scaling  factor  (restoration  stepsize)  in  the 
range  0 < a < 1 . 

Descent  Property.  The  linearized  equations  (33)-(35) 
admit  an  infinite  number  of  solutions.  Each  of  these  solutions 
is  characterized  by  a descent  property  in  the  constraint  error 
P.  This  can  be  seen  by  computing  the  first  variation  of  the 
functional  ( 24 ) : 

5P  = 2 1 (x  - (f)  (A  X-  ())3ax  - ({I'^Au  - (})'^Ati)  dt 
••  0 

rl 

+ 2!  s'^  (S^Ax  + S^Au  + S^Att)  dt  + 2 [ ( '^^  Ax  + i/^^Att  ) ]^  (36) 

Jo  ' X IT 

and  by  observing  that,  when  the  perturbations  defined  by  (33)-(35) 
are  employed,  the  first  variation  of  the  constraint  error  (36) 
becomes 

6P  = -2aP.  (37) 

Since  P > 0,  Eq.  (37)  shov>7S  that  6P  < 0.  Hence,  for  a sufficien- 
tly small,  a decrease  in  the  constraint  error  P is  guaranteed. 
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Auxiliary  Minimization  Problem.  Since  Eqs.  (33)-(35)  ad- 
mit an  infinite  number  of  solutions,  an  additional  requirement  must 
be  introduced  in  order  to  uniquely  define  the  restorative  itera- 
tion. More  specifically,  we  require  that  restoration  be 
accomplished  with  the  least-square  change  of  the  control  and 
the  parameter  (Ref.  10).  Thus,  we  seek  the  minimum  of  the 
quadratic  functional 


K = 


(l/2a) 


C T 

I Au  Audt  + 
^ 0 


T 

An  Ait 


(38) 


with  respect  to  the  perturbations  Ax(t),  Au(t),  An  which  satisfy 
the  linearized  constraints  (33) -(35). 

Special  Variations.  From  calculus  of  variations,  we  know 
that  the  problem  represented  by  ( 33 ),  ( 34),  ( 3 5) , ( 38 ) i s one  of  the  Bolza 
type.  In  this  connection,  lot  A(t)  denote  a variable  Lagrange 
multiplier  associated  with  the  differential  constraint  (33) ; 
let  f (t)  denote  a variable  Lagrange  multiplier  associated  with 
the  nondifferential  constraint  (34);  and  let  y denote  a constant 
Lagrange  multiplier  associated  with  the  final  condition  (35-2) . 
With  this  understanding,  the  Euler  equations  optimizing  Ax(t), 
Au(t),  Ati  and  the  natural  condition  arising  from  the  transver- 
sality  condition  are  written  as 


X + 


(ti  A — 

X 


S p 
X 


= 0, 


0 < t < 1, 


(39) 
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Au=a(4i  X-S  p), 

^u  u 

0 < t < 1, 

(40) 

r 

1 

Att  = a 

l (<^^A  - S^p)dt  - 

'-•'o 

(41) 

(A  + 

o 

tl 

(42) 

Summarizing,  we  seek  variations  Ax(t),Au(t),  Arr  and  multipliers 
A(t),  p(t),  p which  satisfy  the  constraints  (33)-(35)  and  the 
optimality  conditions  (39)-(42). 

Basic  Functions.  If  the  definitions  (29)  are  invoked, 
the  stepsize  a can  be  eliminated,  and  Eqs , (33) -(35)  and 

(39)-(42)  can  be  rewritten  as 


A - 4'^A-  (X-  = 0, 

0 < 

t ■_ 

1, 

(43) 

T T T 

S A + S B + S C + S = 0, 

X u n 

0 < 

t < 

1, 

(44) 

(A)  Q = 0,  (i!^A  + v'J'c  + (i)^  = 0 , 

(45) 

A + ^j^A  - S^p  = 0, 

0 < 

t < 

1, 

(46) 

B - . S^p  = 0, 

0 £ 

t < 

1, 

(47) 

fl 

c + + s^r)dt  + 1 = 0, 

(48) 

•’  0 

= 0- 

(49) 

Equations  (43)-(49)  uniquely  define  the  basic  functions  A(t), 
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B(t),  C as  well  as  the  multipliers  \{t),  p(t),  u of  the  restora- 
tion phase. 

Solution  Technique.  Let  y denote  the  (n+p) -vector 


y = 


A (0) 

c 


(50) 


Let  a sweep  be  defined  as  a forward  integration  of  the  system 
(43) -(49)  obtained  by  (a)  assigning  a particular  value  to  the 
vector  y,  (b)  employing  Eqs.  (43),  (44),  (45-1),  (46),  (47),  and 
(c)  bypassing  Eqs.  (45-2),  (48),  (49). 

Let  n + p + 1 independent  sweeps  be  executed.  More  speci- 
fically, the  first  n + p sweeps  are  executed  by  choosing  the 
n + p linearly  independent  vectors  y to  be  the  columns  of  the 
identity  matrix  of  order  n+  p.  The  last  sweep  is  executed  by 
choosing  y as  the  null  vector.  In  this  way,  we  obtain  the 
particular  solutions  (Refs.  11-14) 


A^  ( t),  (t),  C^,  (t),  (t),  i = 1,2,  ...,n  + p + 1.  (51) 

Then,  we  introduce  the  n + p+1  undetermined,  scalar  con- 
stant and  form  the  linear  combinations 


A (t)  = E)t . A^  (t). 


B(t)  = ?;)<.  B.  (t). 


C = Ik.C.  , 
1 1 


(52) 


A (t)  = Ek.  A , (t). 


p (t)  = (t). 


(53) 
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where  the  summation  is  taken  over  the  index  i.  The  n + p+  1 
coefficients  and  the  q components  of  the  vector  y are 

obtained  by  forcing  the  linear  combinations  (52) -(53)  to 
satisfy  Eqs.  (45-2),  (48),  (49),  together  with  the  normalization 
condition  (Refs.  11-14) 


T.k^  = 1.  (54) 

Stepsize . With  the  basic  functions  A(t),  B(t),  C known, 
we  consider  the  one-parameter  family  of  solutions  (30) . For 
this  one-parameter  family,  the  constraint  error  (24)  becomes  a 
function  of  the  form 


P = P(a)  . (55) 


Then,  the  stepsize  cx  must  be  selected  so  that  the  inequality 

P(ct)  < P(0)  (56) 


is  satisfied  while  keeping 


1 (a) > 0 


(57) 


Satisfaction  of  Ineqs.  (56)  and  (57)  is  guaranteed  for  a suffi- 
ciently small.  Any  violation  of  the  above  inequalities 
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necessitates  a reduction  in  the  stepsize.  Such  a reduction  can 
be  obtained  by  eniployiny  a bisection  process,  starting  from  the 
reference  stepsize 


«Q  = 1 . (58) 

This  reference  stepsize  has  the  following  property;  it  yields 
one-step  restoration  for  the  limiting  case  where  the  constraints 


(13)-(15) 


are  linear. 
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5 . Conjugate  Gradient  Phase:  General  Case 

As  stated  before,  the  conjugate  gradient  phase  is  started 
only  when  the  constraint  error  P satisfies  Ineq.  (27-1).  It 
involves  a single  iteration,  which  is  designed  to  decrease  the 
value  of  the  functional  I or  the  augmented  functional  J,  while 
satisfying  the  constraints  to  first  order.  During  this  itera- 
tion, the  first  variation  of  the  functional  I is  minimized, 
subject  to  the  linearized  constraints.  The  minimization  is 
performed  over  the  class  of  variations  of  the  control  and  the 
parameter  which  are  equidistant  from  some  constant  multiple  of 
the  corresponding  variations  of  the  previous  conjugate  gradient 
phase . 

For  the  sake  of  clarity,  the  general  structure  of  the  con- 
jugate gradient  phase  is  given  f irst  in  this  section . The  linear 

case  is  treated  in  Section  6,  and  the  linear-quadratic  case 
is  treated  in  Section  7.  Then,  the  extension  to  the  nonlinear- 
nonquadratic  case  is  given  in  Section  8. 


Linearized  Equations.  Let  x(t),  u(t),  it  denote  nominal 
functions  satisfying  ( 13  ) - (1 5),  ^ and  let  x(t),  u(t),  tt  denote 


^These  nominal  functions  can  be  obtained  by  employing  the  res- 
toration algorithm  of  Section  4. 


27 


AAR-126 


varied  functions  also  satisfying  (13) -(15).  To  first  order, 
the  perturbations  Ax(t),  Au(t),  An  must  satisfy  the  linearized 
constraint  equations 


• 'p  rp  (p 

Ax  - (Ji  Ax  - (()  Au  - 4)  An  = 0, 
X u n ' 


O',  t < 1, 


T T T 

S Ax  + S Au  + S A7r  = 0, 
X u n 


0 < t •:  1, 


(Ax) Q = 0, 


( + i/j'^An  ) = 0 . 

X ^n  1 


Note  the  difference  between  Eqs.  (33) -(35)  and  Eqs.  (59) -(61). 
While  the  former  are  nonhomogeneous , the  latter  are  homogeneous 
in  the  perturbations  Ax(t),  Au(t),  An. 

Auxiliary  Minimization  Problem.  Since  the  linearized 
equations  (59) -(61)  admit  an  infinite  number  of  solutions, 
some  additional  requirement  must  be  introduced  in  order  to 
uniquely  define  the  conjugate  gradient  iteration.  More  spe- 
cifically, we  consider  the  first  variation  of  the  functional 
(12)  : 


<Sl=\  (f  Ax+f  Au  + f^.  .■n)dt  + (g^Ax  + g^An  )t 
I X u n ^x  -’n  l 

J n 


and  the  isoperimotr ic  constraint 


K = I ( Au  - (=  Au  ) "^  ( Au  - (’i'  u)  rit  + ( Ar  - BAn  ) "^  ( An  - B Att  ),  ( 63 ) 


r 
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1 


where  K and  6 are  undetermined  constants.  The  symbols 
Ax(t),  Au(t),  Att  denote  the  variations  associated  with  the 
previous  conjugate  gradient  iteration.  Therefore,  because 
of  ( 29 ) , we  have 

Ax(t)=aA(t),  Au(t)=aB(t),  ATT=aC.  (64) 

Then,  we  seek  the  minimum  of  the  linear  functional  (62)  with 
respect  to  the  perturbations  Ax(t),  Au(t),  Air  which  satisfy  the 
linearized  constraints  (59) -(61)  and  the  quadratic  isoperi- 

metric  constraint  (63) . i 

Special  Variations.  From  calculus  of  variations,  we 
know  that  the  problem  represented  by  (59) -(63)  is  one  of  the 
Bolza  type  with  an  added  isoperimetric  constraint  on  the 

1 

variations  of  the  control  and  the  parameter.  In  this  connec-  ; 

tion,  let  X(t)  denote  a variable  Lagrange  multiplier  associated 
with  the  differential  constraint  (59);  let  p(t)  denote  a 

variable  Lagrange  multiplier  associated  with  the  nondifferen-  ! 

I 

tial  constraint  (60);  let  m denote  a constant  Lagrange  1 

multiplier  associated  with  the  final  condition  (61-2);  and  let  , 

l/2a  denote  a constant  Lagrange  multiplier  associated  with  the 

isoperimetric  constraint  (63).  With  this  understanding,  the 

Euler  equations  optimizing  Ax(t),  Au(t),  Ati  and  the  natural 

condition  arising  from  the  transversal ity  condition  are 


written  as 
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A - f +4'  A - s 

X X X 

p = 0 , 

0 < t < 1 , 

(65) 

(Au  - BAu)/a  + 

f,,  - 4,  A + S p = 0, 

0 < t < 1, 

(66  ) 

u u u 

~ — 

1 

• 1 

(Att  - PAtt)/ a + 

(f  -4  A + S p)dt  + 

7T  7T  7T 

(67  ) 

0 

(A  + g^+4^u)l  = 

0 . 

(68  ) 

Summarizing,  we  seek  variations  Ax(t),  Au{t),  An  and  multipliers 
A(t),  p(t),  p,  l/2u  which  satisfy  the  constraints  (59),  (60),  (61), 
(63)  and  the  optimality  conditions  (65)-(68). 

Basic  Functions.  Let  the  definitions  (29)  be  invoked 
for  both  the  present  conjugate  gradient  phase  and  the  previous 
conjugate  gradient  phase.  Let  the  directional  coefficient y 
be  defined  as 


Y = B (a/a)  . (69) 

With  this  understanding,  the  stepsize  a can  be  eliminated,  and 
Eqs.  (59)-(61)  and  (65)-(68)  can  be  rewritten  as 

A - ^,^A  - (()^B  - <))^C  = 0,  0<t<l,  (70) 

S^’a+  S^B  + S^C  = 0,  0<t<l,  (71) 

(ip'^A  + iJ/^C)^  = 0, 


(A)q  = 0, 


(72) 
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X-f  +(})X-Sp=0, 

X X X 

B — '("B+f  — X + S p=0, 

u u u 

C-YC+  \ (f  -4)  X + S p)dt 
^ 0 

(X  + = 0 . 


0 < 

t < 1, 

(73) 

0 < 

t < 1, 

(74) 

(g  + = 

^ TT  TI  1 

0, 

(75) 

(76) 

For  a given  value  of  the  directional  coefficient  y , Eqs. 
(70)-{76)  uniquely  define  the  basic  functions  A(t),  B(t),  C as 
well  as  the  multipliers  X(t),  p{t),  p of  the  conjugate  gradient 
phase. 

Isoperimetr ic  Constant.  In  the  light  of  (29)  and (69),  the  iso 
perimetric  functional  (63)  tal<es  the  form 


K = 


(B  - yB)'^(B  - YB)dt  + (C  - yC)'^(C  - yC) 

0 


(77) 


If  the  basic  functions  A(t),  B(t),  C are  consistent  with  (70)-(76) 
the  error  in  the  optimality  conditions  (25)  reduces  to 

Q=\  (B  - yB)'^(B  - YB)dt  + (C  - yC)'^(C  - yC)  . (78) 

Jo 


Consequently,  the  following  relation  ties  the  isoperimetr ic 
constant,  the  stepsize,  and  the  error  in  the  optimality  con- 


ditions : 


K = a^Q  . 


(79) 
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Clearly,  to  assign  values  to  the  isoper imetr ic  constant 
is  the  same  as  assigning  values  to  the  stepsize.  However, 
there  is  no  clear-cut  way  of  determining  a priori  convenient 
values  for  the  isoperimetr ic  constant.  Therefore,  the  imple- 
mentation of  the  conjugate  gradient  algorithm  becomes  simpler 
if  one  avoids  evaluating  a in  terms  of  K through  (79)  and 
assigns  values  to  u directly. 

Descent  Property.  Next,  consider  the  augmented  functio- 
nal (16)  and  its  first  variation 

6J=I  (f  -4>.\  + Sp-A)  "^Axdt  + 1 ( ~ 4',,  A + S p ) "^Audt 

J Q Jo 


1 

rrl 

T 

r T ' 

+ j 

1 \ S^0)dt+  (g^+  ^'^P)p 

Att  + 

L 

When  the  perturbations  defined  by  (29)  and  (70) -(76)  are 
employed,  Eq.  (80)  becomes 


6J  = -a 


c ^ ip 

\ (B  - , B)  Bdt  + 
•^0 


(C  - C)'^C 


(81) 


Upon  invoking  Eq.  (78)  and  defining  the  quantity 

Z=  1 (B-  yB)  Bdt+  (C  - yC)  C , 

0 


(82) 


we  see  that  E(} . (81)  can  ho  rewritten  as 
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6J  = -a (Q  + yZ)  . (83) 

For  the  first  iteration  of  the  conjugate  gradient  phase, 
one  sets 


Y = 0, 


(84) 


with  the  implication  that 


6J  = -aQ  . (85) 

Since  Q > 0,  Eq.  (85)  shows  that  6J<0.  Hence,  for  a sufficien- 
tly small,  a decrease  in  the  augmented  functional  J is 
guaranteed . 

For  subsequent  iterations,  one  sets  , ^ 0.  More  specifi- 
cally, the  directional  coefficient  must  be  such  that 

t ''  0,  (86) 

and  its  proper  value  is  discussed  in  Section  7 . At  any  rate, 
Eq.  (83)  shows  that  6J<0  providing 

Q + YZ  > 0 , (87) 


where  Q is  given  by  (78)  and  Z is  given  by  (82).  Hence,  for 
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3 3 


a sufficiently  small,  the  decrease  in  the  augmented  functional 
J is  guaranteed  as  long  as  Ineq.  (87)  is  satisfied.  If  Ineq. 
(87)  is  violated,  the  descent  property  on  J no  longer  holds, 
and  the  conjugate  gradient  phase  must  be  restarted  by  resetting 
the  directional  coefficient  y at  the  level  (84) . 

Solution  Technique.  Now,  assume  that  a particular  value 
is  given  to  the  directional  coefficient  y.  Let  y denote  the 
(n  + p) -vector 


y = 


A (0) 
C 


(88) 


Let  a sweep  be  defined  as  a forward  integration  of  the  system 
(70) -(76)  obtained  by  (a)  assigning  a particular  value  to  the 
vector  y,  (b)  employing  Eqs.  (70),  (71),  (72-1),  (73),  (74),  and 
(c)  bypassing  Eqs.  (72-2),  (75),  (76). 

Let  n + p + 1 independent  sweeps  be  executed.  More  speci- 
fically, the  first  n+p  sweeps  are  executed  by  choosing  the 
n+p  linearly  independent  vectors  y to  be  the  columns  of  the 
identity  matrix  of  order  n + p . The  last  sweep  is  executed  by 
choosing  y as  the  null  vector.  In  this  way,  we  obtain  the 
particular  solutions  (Refs.  11-14) 


A.  (t),  D.  (t),  C. 


(t),  (t). 


i = l,2,...,n  + p+  l. 


(89) 


w 
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Then,  we  introduce  the  n + p+  1 undetermined,  scalar  cons- 
tants and  form  the  linear  combinations 

A(t)  = Zk^A^  (t),  B(t)  = Zk^B^(t),  C=Zk.C^,  (90) 

A (t)  = Zk^A^  (t),  p (t)  = Zk.p^ (t),  (91) 

where  the  summation  is  taken  over  the  index  i.  The  n + p+1 
coefficients  k^  and  the  q components  of  the  vector  p are 
obtained  by  forcing  the  linear  combinations  (90) -(91)  to  satisfy 
Eqs.  (72-2),  (75),  (76),  together  with  the  normalization  condi- 

tion (Refs.  11-14) 

i:k^  = 1 . (92) 

General  Solution.  Next,  assume  that  two  particular  va- 
lues are  given  to  the  directional  coefficient  7,  for  instance, 

VO  and  Y**  = l.  (93) 

Assume  that  the  previous  solution  technique  is  employed  twice, 
and  denote  by 


A *(t),  B^  (t),  , 


\ (t),  p^  (t),  u* 


(94) 


and 
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A**(t),  D*,(t),  C**,A*^(t),  P**(t),  ;.**  (95) 

the  particular  solutions  of  (70) -(76)  corresponding  to  (93-1) 
and  (93-2),  respectively.  Simple  manipulations,  omitted  for 
the  sal<e  of  brevity,  show  that  the  general  solution  of 
(70) -(76),  valid  for  any  value  of  the  directional  coefficient 
y,  can  be  written  as 

A(t)  = A*  (t)  + Y [A**  (t)  - A*  (t)],  B(t)  = B*(t)  + Y [B**(t)  -B*(t)], 

C = C*  + ',  (C**- C*),  (96) 

A(t)  = A*  (t)  + > [A**  (t)  - A*  (t)],  p(t)  = e*(t)  + ■,  [r**(t)  -p*(t)], 

(J  = + Y (m**-  P*)  . (97) 


As  a conclusion,  the  general  solution  of  (70) -(76)  requires 
that  two  sweeps  of  n + p + 1 integrations  be  executed,  one  lead- 
ing to  the  particular  solution  (94)  and  one  leading  to  the 
particular  solution  (95).  However, if  the  constraints  are  li- 
near, the  general  solution  of  (70) -(76)  requires  only  one  sweep 
of  n+p+1  integrations,  that  leading  to  the  particular  solu- 
tion (94),  as  is  shown  in  Section  6. 


Stepsizc  and  Directional  Coefficient.  After  the  general 
solution  (96) -(97)  is  known,  the  next  step  is  to  determine  the 
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proper  values  of  the  stepsize  a and  the  directional  coefficient 
> . A logical  scheme  is  that  of  determining  these  quantities 
so  that  the  augmented  functional  (16)  is  minimized. 

For  the  varied  functions  x{t),  u(t),  tt  , let  the  augmented 
functional  (16)  be  written  in  the  form 

- T~  ^ T'  T-  -T'  - 

J=  (X  x)  +\  (f  - X 0 + P S -\  x)dt+  (g  + p •^)  . (98) 

u J 0 

In  view  of  (30)  and  (96),  the  varied  functions  x(t),  u(t),  ti  cons- 
titute the  two-parameter  family 


X(t) 

= x(t)  + a ^A*  (t)  + Y [A^,*  (t)  - A^(t)l  1 

1 

(99-1) 

u(t) 

= u (t)  + a (t)  + Y [B**  (t)  -B*  (t)  ] I 

1 

(99-2) 

Tt 

= 7r  + a[C^+y(Cj^^-C^)]  . 

(99-3) 

On  the  other  hand,  the  multipliers  A(t),  p(t),  p constitute  the 
one-parameter  family  (97).  Upon  using  (97)  and  (99),  we  see 
that  the  augmented  functional  (98)  takes  the  form 

J = J(a,^)  . (100) 

Therefore,  the  optimum  values  of  a and  y satisfy  the  relations 


r 
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( a , Y ) = 0 , 


(a,Y)  = 0 . 


(101) 


In  principle,  one  can  solve  Eqs.  (101)  using  an  exact 
search,  as  done  in  Ref.  15  for  mathematical  programming  prob- 
lems. The  resulting  algorithm  constitutes  the  extension  to 
optimal  control  problems  of  the  memory-gradient  method  of 
Ref.  15.  However,  the  simultaneous  determination  of  a and  y might 
be  expensive  computationally;  this  being  the  case,  we  follow 
a different  road.  First,  we  determine  an  approximate  value  of 
the  directional  coefficient  , based  on  the  consideration  of 
the  linear-quadratic  case  (Section  7).  Once  *,  is  )cnown,  the 
two-parameter  family  (100)  reduces  to  the  one-parameter  family 


J = J (a) 


Then,  the  optimum  stepsize  n satisfies  the  relation 


J^(a)  = 0, 


(102) 


(103) 


whose  numerical  solution  can  be  arrived  at  in  a variety  of 
ways.  For  example,  within  the  frame  of  the  linear-quadratic 
case,  the  numerical  solution  of  (103)  can  be  obtained  with 
quadratic  interpolation  (Section  8 ) . On  the  other  hand,  within 
theframeof  the  nonlinear-nonquadratic  case,  the  numerical 


I 


38 


AAR-126 


solution  of  (103)  can  be  obtained  with  cubic  interpolation 
(Section  3 ) . 


39 


aar-126 


6.  Conjugate  Gradient  Phase:  Linear  Constraints 

In  the  previous  section,  we  derived  some  general  rela- 
tions which  are  valid  for  the  conjugate  gradient  phase, 
regardless  of  the  analytical  form  of  the  functional  (12)  and 
the  constraints  (13) -(15).  In  this  section,  we  give  the  parti- 
cular relations  whicli  are  valid  if  the  constraints  arc  linear. 

General  Solution.  Under  the  linearity  assumption  for  the 
constraints,  consider  the  system  (70) -(76)  which  defines  the 
basic  functions  A(t),  B(t),  C as  well  as  the  multipliers  A(t), 
p(t),  u.  By  substitution,  it  can  be  verified  that  the  particu- 
lar solutions  (94)  and  (95)  satisfy  the  relations 


^**(t)-A*(t)  = A(t), 

(t)-B^(t)  =B(t), 

c**- 

- c = c 

(104) 

>**(t)-  \^(t)=  0, 

P**(t)-  P*(t)  = 0, 

P** 

1 

II 

o 

(105) 

As  a consequence,  Eqs, 

. (96) -(97)  take  the 

simpler 

form 

A(t)  = A*(t)  + YA(t), 

B(t)=  B*  (t)  +-rB(t), 

+ 

■K 

u 

II 

u 

iC, 

(106) 

Mt)  = (t). 

p(t ) = p*  (t). 

11=0* 

. 

(107) 

Tlu'  imjilication  of  (106) -(107)  is  the  following:  under  the 
assumption  of  liiuar  consti'a  i nt  s , the  general  solution  of 
(7  0) -(7())  (Mn  now  be  obtained  by  executing  only  one  sweep 
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(instead  of  two)  of  n + p + 1 integrations,  namely,  the  sweep 
necessary  to  generate  the  particular  solution  (94)  . This  is 
the  solution  corresponding  to  (93-1),  namely,  the  solution 
associated  with  the  ordinary  gradient  method  of  Ref.  9. 

Isoperimetric  Constant.  Under  the  linearity  assumption 
for  the  constraints,  Eq. (79)  still  holds,  but  the  error  in 
the  optimality  conditions  (78)  simplifies  to 

0 = I B*B*dt  + cjc*  . (108) 

J 0 

Descent  Property.  Under  the  linearity  assumption  for 
the  constraints,  Eq. (83)  still  holds,  but  the  functional  (82) 
simplifies  to 

Z = I B*Bdt  + C*C  . (109) 

0 


Local  Orthogonality  Conditions.  Under  the  linear- 
ity assumption  for  the  constraints , the  two-parameter  family  (99) 
simplifies  to 


x(t)  = x(t)  + a[A*(t)  + YA(t)], 
u(t)  = u(t)  +fx[B,,(t)  +YB(t)l, 


(110-1) 

(110-2) 


71 


71  -*•  a (C*  + yO  . 


(110-3) 


r 
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Next,  we  consider  the  augmented  functional  (98)  and 
observe  that,  for  the  two-parameter  family  (110),  it  takes  the 
form  (100).  Therefore,  the  optimal  values  of  the  stepsize  a 
and  the  directional  coefficient  y satisfy  the  relations  (101), 
Because  of  the  assumed  linearity  of  the  constraints  (13)-(15), 
Eqs.  (101)  take  the  particular  form 


L 


iV*. 

J 0 ^ 


A + S p - A ) '^Ad  t + I ( f - A + S p ) ^ Bd  t 
X I u ' u u 


’ c 

T 

r-  - 

T 

-I- 

L 0 J 

C + 

=0,  (111) 


f (f  ’ 

I X X 

J 0 


A+S  p-A)'^Adt+l  (f  -<}>  \ + S (-))'^Bdt 
X ) u u u 

J 0 


' r ^ ~ 

T 

- 

1 (f  —0  A+S  p)dt  + (g  + 0 u)i 

1 TT^iT  71  -’it  ^n’  l 

C + 

{ A + kj  + ;i ) A 

-’x  X 

L-^0  J 

=0,  (112) 


with  the  implication  that 


^ T 

(J'  A + S p - i ) A*dt  + 1 (f  - \ + s p)  '^B*dt 

„ X ^ X X I u ' u u * 

0 0 


Cl 

T 

^ T 

\ +S^P)dt+(g^  + 0^ti)j 

L •'0 

c*  + 

(A+ g^+  0^11)  A* 

=0.  (113) 


Bee  luse  of  tlie  linearity  of  the  constraints  (13) -(15)  and  after 
invoking  Eqs.  (70) -(72),  one  can  show  that  Eqs.  (Ill) -(1)3) 
hold  for  any  distribution  of  Lagrange  multipliers.  In  narticu- 
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lar,  they  hold  if  (t),  p(t),  ',i  arc  replaced  by  A(t),  p(t),  u- 
This  yields  the  supplementary  relations 


j ( - X Adt  + I 


p-X)’^Adt+(  + S^P^Bdt 

0 


- 

T 

■ . . . T ' 

\ (f  -ifi  A + S p)  dt  + 

1 TT  ^7T 

L-’O 

^TT+  ^^^l 

c 

(X+g^+ii'^u)  A 

= 0,  (114) 


C' 

[s: 


{ f - ({!  \ + S p 
X x^ 


+ \ (f„-  ; \+  S p)dt 

‘ 1 TT  ^ V TT 

'o 


!}’■- 

[i 


-+  \ A+S  p)dt+  (c7  + il' 

It  tt’tt  it  -bi^ 


' T r 

- \ ) Adt  + \ 

+ { g + P ) T 
TT  ^ TT  1 

■'  0 


“t 


4'  X + S o - A ) A*dt  + 

^ X 


(X  + g^+  A 


+ S^n)  B*dt 


=0,  (115) 


(X  + g^+  i'^L)  A* 


= 0.  (116) 


Next,  we  combine  Eqs.  (114) -(116)  with  Eqs.  (73) -(76)  written 
for  the  next  iteration.  This  leads  to  the  following  local 
orthogonality  conditions: 


B^Bdt  + C*C  = 0, 

-w  w AS  ^ rn  •. 

B*Bdt  + C^C  = 0, 


(117-1) 


(117-2) 


(117-3) 
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Here,  the  adjective  local  is  employed  to  mean  that  Eqs.  (117) 
involve  vectors  B{t),  C which  are  solutions  of  (70) -(76)  compu- 
ted for  the  present  iteration  and  the  previous  iteration;  they 
also  involve  vectors  B^(t),  C*  which  are  solutions  of  (70) -(76) 

for  = 0 computed  for  the  present  iteration  and  the  next  iter- 
ation. 
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4 A 


7 . Conjugate  Gradient  Phase;  Quadratic  Functional  and 

Linear  Constraints 

In  the  previous  section,  we  assumed  that  the  constraints 
(13)- (15)  are  linear  and  arrived  at  the  local  orthogonality 
conditions  (117).  in  this  section,  we  retain  the  constraint 
linearity  hypothesis  and  further  asssume  that  the  functional 
(12)  is  quadratic. 

For  the  sake  of  compactness,  let  y and  E denote  the 
vectors 


(118) 


Let  f^  and  g^  denote  the  gradients  of  the  functions  f and  g 
with  respect  to  the  vector  y; 


(119) 


Under  the  assumption  that  the  functions  f and  g are  quadratic 
in  their  respective  arguments,  the  following  exact  relations 
can  be  established: 


f 

y 


+ (if  E , 

yy 


g + aq  E 

y yy 


(i?.o) 
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wher 


” 

“ 

~ 

- 

f 

XX 

f 

XU 

c 

■^XTT 

g 

^xx 

0 

g 

X7I 

l-tl 

II 

f 

ux 

f 

uu 

f 

UTT 

U3 

II 

0 

0 

0 

f 

nx 

f 

••U 

f 

TTTT 

g 

^TTX 

0 

97T7T_ 

With  this  understanding,  Eqs.  (Ill) -(113)  become 

1 ,.1 


■ ^)"Adt+  1 (f  - .j;  A+S  p)'^’Bdt 

0 X X X a u u 


IS 


1 


+ 1 \ (f,-  + S^p)dt+  (g^+ 

’o 


C + 


(\+  g^+  ^A 


+a  \ \ E'^f  Edt  + (b'^q  E),  = 0, 

I )q  yy  - yy  i 


1 


1 . fi 


^^x“  '^'x'  - A)  ^Adt  + (fu" 

0 •'0 


+ \ (f  -<{■  A+S  p)dt+(g  + !'  u), 

I*  n^TT  71  -"ti  7T1 


[S. 

■is: 


'0 

•1 


^ /N 

C + 

(A+g^+ '-Jj^m)  a 

ta  \ E^f  Edt  + (E^g  E) , = 0, 

II  yy  Jyy  ' 1 

'0  -I 


S. 

[S 

■[s:"''v7^ 


^^x"  '^x'^  ■*■  ^xP  “ A)'^A*dt+  \ (f^-  <i^A  + S^p)'^B*dt 

•'0 


'^r  •7i'^+  )dt+  (cj^+>i7^ii)^ 


c*+ 


( '+  gj^+ 


:*dt+  (E 


(121) 


1 

(122) 


(123) 


0. 


(124) 
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Upon  invoking  Eqs.  (73)-(76),  we  see  that  Eqs.  (122)-{124)  can 
be  rewritten  as 


T TIT  T 

B*Bdt  + Cic-cx  \ E f Edt+  (E  q E), 

* * LJq  yy  ^yy  i 


(, 

B^Bdt  + C^C-aj^J 


Edt  +(E  g E), 
' ^vy  1 


= 0, 


= 0, 


dt  + (E  g E*)  , 
^yy  * 1 


= 0, 


(125-1) 


(125-2) 


(125-3) 


Local  Conjugacy  Condition.  Next,  we  employ  the  local 
orthogonality  condition  (117-1)  written  for  the  present  itera- 
tion, and  observe  that  (125-2)  yields  the  local  conjugacy 
condition 


I £”^1  Edt+(E'^g  E),  =0.  (126) 

I yy  ^yy  '1  V J-  u / 

0 

Here,  the  adjective  local  is  employed  to  mean  that  Eq.  (126) 
involves  vectors  E(t),  that  is,  vectors  A(t),  B(t),  C,  which  are 
solutions  of  (70)- (76)  computed  for  the  present  iteration  and 
the  previous  iteration. 

Stepsize.  After  observing  that 

E=E*+,E  (127) 

and  making  use  of  the  local  conjugacy  condition  (126),  we  see 
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that  Eq.  (125-3)  can  be  rewritten  as 


BlB*c]t  + I \ Edt  + (E^g,..  E), 


T 

" C 


yy 


^yy  '1 


= 0, 


(128) 


This  equation  enables  one  to  compute  the  optimum  stepsize  a, 
once  the  value  of  the  directional  coefficient  y is  known. 

Directional  Coefficient.  After  involving  Eq . (127),  we  see 

that  the  local  conjugacy  condition  (126)  becomes 


1 


B*f,,.,Edt+  (Eig„,  E),  + V 


yy 


'yy  '1 


rr-- 

\ E^f  Edi 

L3o 


" ' ‘ """  Edt  + (E^g  E),  , 

"^yy  ij 


0.  (129) 


'0  '-•^0 
If  we  employ  Eq.  (128)  written  for  the  previous  iteration, 
Eq.  (129)  becomes 

rd 


' ^ip  ^ I I ~ T ' T 

B*D*dt  + +a|  \ E*f„,,Edt+  (E*g,„,E). 


yy 


yy  '1 


= 0 


'0  J L--  0 

and,  in  the  light  of  Eqs.  (120),  can  be  rewritten  as 


y I B^B*dt  + C 


= 0, 


(130) 


(131) 


If  Eqs.  (70)-(76)  and  (117)-(119)  arc  employed,  the  follow- 
ing relations  can  be  shown  to  hold: 


1' 


dt  + 


y 


( bIb. 

j 0 


dt+  C*CA  = 0 


(132) 
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E:f^dt+  (E*g^)^  = 0 . 


(133) 


As  a consequence,  Eq.{131)  becomes 


^ I B^B*dt  + c:c* 


BlB*dt+cIc*  = 0 


(134) 


and  can  be  rewritten  as 


= Q/Q  , 


(135) 


where 


1 rp  ip 

Q=  \ B*B*dt  + C*C*  , 


j 0 

fl.T 

•1  B. 

J n 


A rri 

B*dt  + C*C* 


(136) 


(137) 


denote  the  errors  in  the  optimality  conditions  for  the  present 
iteration  and  the  previous  iteration,  respectively.  These 
quantities  are  known,  since  they  involve  vectors  B^(t),  C* 
which  are  solutions  of  (70)-(76)  for  y = 0 computed  for  the 
present  iteration  and  the  previous  iteration. 

Descent  Property.  Because  of  the  local  orthogonality 
condition  (117-1)  written  for  the  previous  iteration,  Eq . (109) 

yields 


Z = 0. 


L 
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As  a consequence,  the  first  variation  of  the  augmented  fun- 
ctional (83)  reduces  to 


'"J  = - tQ  , (139) 

whore  the  error  in  the  optimality  condition  Q is  given  by  Eq. 
(136).  Equation  (139)  holds  for  any  conjugate  gradient  iter- 
ation and  shows  that,  since  Q • 0,  we  have  ■ 0.  Hence,  for 
n suffic  ently  small,  a decrease  in  the  augmented  functional  J 
is  guaranteed.  Tn  conclusion,  for  the  linear-quadratic  case, 
the  restart  procedure  mentioned  in  Eection  5 never  occurs. 

This  means  that  the  directional  coefficient  is  set  at  the 
level  (84)  only  for  tlie  first  iteration. 

Genera 1 Or thogonality  and  Conjugacy  Cond i t io ns . Assume 
now  that  the  algorithm  described  by  Eqs.  (70) -(76)  and  (110) 
is  employed,  starting  with  a feasible  nominal  solution.  Fur- 
ther, assume  that  the  first  conjugate  gradient  iteration  is 
done  with 


Y = 0,  (140) 

meaning  that  this  is  an  ordinary  gradient  iteration.  Under 
these  assumptions  and  with  reference  to  the  linear-quadratic 
case,  one  can  generalize  the  local  orthogonality  conditions 
(117)  and  the  local  conjugacy  condition  (126)  as  follows: 
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and 


BlBpdt  + cJCp=  0 , 


'0 


B*B*pdt  + C^C^p  0 , 


1 


E'^f  E dt  + (e’^q  E ),  = 0 , 

yy  p ^yy  p'l 


(141-1) 

(141-2) 

(142) 


where  the  subscript  p denotes  any  iteration  preceding  the 
present  iteration.  While  these  equations  do  not  guarantee 
convergence  in  a finite  number  of  steps,  they  do  guarantee  that 
the  algorithm  generates  a sequence  of  linearly  indepentent 
vectors  E(t),  that  is,  a sequence  of  linearly  independent  va- 
riations per  unit  stepsize  A(t),  B(t),  C. 


5 
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8.  Conjugate  Gradient  Phase:  Practical  Implementation 

In  this  section,  we  summarize  the  results  of  Sections 
5-7,  and  suggest  practical  ways  of  utilizing  these  results, 
while  avoiding  tho  use  of  second  derivatives.  This  is  an 
essential  characteristic  of  a first-order  method. 

-Auxiliary  Functions.  The  first  step  is  to  solve  Eqs. 
(70) -(76)  for  a fictitious  value  of  the  directional  coeffi- 
cient, namely. 


. * = 0 . (143) 

This  yields  the  following  linear,  two-point  boundary-value 
problem: 


. T T T 

A.  - ^ C.  = 0, 

0 < t < 

1, 

(144) 

T T T 

S A.  + S B.  + S C.  = 0, 

0 < t < 

1, 

(145) 

(AJ^=C, 

Tc*>l  = “■ 

(146) 

'^x^*  " ^x*^*  ""  ° ' 

0 < t < 

1, 

(147) 

B*+f  -if  + S P*=0, 

* u u * u * 

0 < t 

1, 

(148) 

' 1 

S^pjdt  + (g  + V U*), 

= 0, 

(149) 

0 

(150) 
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Using  the  solution  technique  of  Section  5,  we  obtain  the  fol- 

O 

lowing  auxiliary  functions  and  multipliers  : 

A*(t),  B*(t),  C*,X*(t),  p*(t),  p*  . (151) 

Directional  Coefficient.  The  second  step  is  to  compute 
the  actual  value  of  the  directional  coefficient 
■'1  . For  the  first  conjugate  gradient  phase,  we  set 

Y = 0,  (152) 

meaning  that  the  conjugate  gradient  iteration  is  an  ordinary 
gradient  iteration.  For  subsequent  conjugate  gradient  phases, 
we  set 

Y = Q/Q  , (153) 

where 

Q=  I B*B*dt  + C*C*  , (154-1) 

J 0 

Q=  I B*B*dt  + C*C*  (154-2) 

J 0 

denote  the  errors  in  the  optimality  conditions  at  the  beginning 


g 

These  functions  and  multipliers  are  identical  with  those 
solving  the  linear,  two-point  boundary-value  problem  asso- 
ciated with  the  ordinary  gradient  phase  of  Ref.  9. 
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of  the  present  conjugate  gradient  phase  and  at  the  beginning 
of  the  previous  conjugate  gradient  phase,  respectively. 

Note  that  the  directional  coefficient  (153)  is  accepta- 
ble only  if 

J ^ (0) < 0 , (155) 

where  J^(0)  is  given  by  hq . (160).  if  Ineq.  (255)  is  viola- 

ted, then  the  directional  coefficient  (153)  i.iust  be  discarded 
and  replacea  by  the  value  (152)  . This  i.cans  that  the  algorithr.i 
must  be  restarted  by  replacing  tli*^  conjugate  gradient  phase 
with  an  ordinary  gradient  phase. 

Basic  Functions.  The  third  step  is  to  compute  the  basic 
functions  A(t),  B(t),  C and  the  multipliers  X (t),  p(t),  )i . This 
is  done  with  the  following  formulas: 


A(t) 

= A*  (t)  + ,A(  t). 

B(t) 

= B*  (t)+  YB(t), 

C = C*+  ,C, 

(156) 

(t) 

= ^.(t), 

r-  (t) 

= .T*  (t). 

U = 11*  . 

(157) 

Stepsize.  With  the  basic  functions  A(t),  B(t),  C known, 
we  consider  the  one-parameter  family  of  solutions  (30).  For 
this  one-parameter  family,  the  augmented  functional  (16)  takes 
the  form 

J(u)=(>'^i)j+  f (f  - + p'^S-  i'^x)dt  + (g  + , 

J 0 


(158) 


54 


AAR-126 


with  the  implication  that 


!,(«)=  (f^- 

J n 


A + S p - A ) Adt  + 

X X 


C-f 

1 u ^ u 
J n 


A + S p)  Belt 
u 


rr-  ~ ~ ~ ~ ' 

I' 

- _ 

c + 

(A  + 5 + ^ U)'^A 

A,  X 

, (159) 


and  with  the  further  implication  that 


J^(0)  = -(Q+  YZ)  , 


:i60) 


where 


I rp  rp 

Q=  \ B*B*dt+C*C*  , 


:i6i) 


■1 


z = \ B;Bdt  + c;c 


(162) 


Note  that 


Z=0  (163) 

in  the  linear-quadratic  case  and  that 

2/0  (164) 

In  the  general  case. 

In  addition  to  the  augmented  functional,  the  constraint 
error  (24)  must  be  monitored  during  the  conjugate  gradient 
phase.  For  the  one-parameter  family  of  solutions  (30),  the 
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functional  (24)  takes  the  form 


P(u)=  ( N(x-^)dt+rN(S)dt 

J n Jo 


+ N (!);)- 


(165) 


where  N(b)  denotes  the  norm  squared  of  the  vector  b. 

With  these  preliminaries  in  mind,  the  stepsize  u must 
be  selected  so  that  the  following  inequalities  are  satisfied; 


J^(  ■)/J^(9) 


- 3 ' 


(166) 


J(  ) J(0)  , 


subject  to 


P (a)  '■  ?*  , 


T (a)  • 0 , 


(168) 


where  and  P*  are  preselected  numbers.  While  is  a 
small  number,  P*  need  not  be  necessarily  small. 

Quadratic  Interpolation.  In  the  linear-quadratic  case, 
satisfaction  of  (166) -(167)  can  be  achieved  by  using  quadratic 
interpolation.  The  procedure  is  as  follows. 

Let  the  function  (158)  be  written  in  the  quadratic  form 


J ( 0 - kg  + + k2a  , 


(169) 


with  the  implication  that 


J^(.)  .2k2.  , 


J (\)  = 2k-  . 

K'i  2 


(170) 
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Therefore,  the  coefficients  and  are  given  by 

kQ  = J(0),  k^  = J^(0)  . (171) 

The  coefficient  k2  can  be  computed  by  evaluating  the  augmen- 
ted functional  at  some  . eferencc  stepsize,  for  instance,  u=l. 

If  this  is  done,  one  concludes  that 

k,=J(l)-J(0)-J{0).  (172) 

z (X 

With  the  coefficients  known,  the  optimum  value  of  u can  be 
computed  from  the  relation 

(a)  = 0 , (173) 

which  implies  that 

ao  = -k^/2k2  . (174) 

In  the  linear  quadratic  case,  the  representation  of  the 
augmented  functional  (158)  by  means  of  the  quadratic  form 
(169)  is  exact.  Therefore,  the  quadratic  interpolation  pro- 
cess is  employed  only  once  (one-step  quadratic  interpolation) . 

Cubic  Interpolation.  In  the  general  case,  it  is  better 
to  try  to  achieve  satisfaction  of  (166) -(167)  with  cubic 
interpolation.  The  procedure  is  as  follows. 

We  consider  the  reference  stepsize  k and  the  sequence 


I 


of  stepsizes 
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{xj  = {O,  k,  2k,  4k,  ak,..,}  . (175) 

For  every  element  of  the  sequence  (175) , we  compute  the  quan- 
tities (158) -(159).  Vi/e  denote  by  and  1X2  the  smallest 

consecutive  elements  in  the  sequence  (175)  such  that  the  fol- 
lowing inequalities  are  satisfied 

0 , (a^)  " 0 . (176) 

Then,  assuming  that  the  derivative  J^_^(u)is  continuous,  a relative 
minimum  of  J(n)  occurs  for  a value  such  that 

1 * n ^ 2 ■ (177) 


In  order  to  find  the  minimum  of  J(u)  numerically,  we  ap- 
proximate the  function  J(a)  with  the  cubic  form 


J (a)  = , 


(178) 


with  the  implication  that 


J^(a)=kj^  + 2k2a+3k2a  , 


J^,(a)=2k.,  + 6k^ot  . 
Ota  2 3 


(179) 


The  coefficients  k^  are  computed  by  forcing  the  cubic  function 
(178)  and  its  derivative  to  satisfy  the  exact  values  of  the 
ordinate  J ( t)  and  the  slope  J,^(a)  at  and  ; that  is,  the 

coefficients  k^  are  computed  from  the  conditions 


L 
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J( 

= kQtk^Ui  + k2a2  + k3.tJ  , 

(180-1) 

J (a2) 

= kg  + ^^1X2  + ^2**2  ^ ' 

(180-2) 

) = k^^  + 2k20t^  + Bk^a^  , 

(180-3) 

) = k^  + 2k2U2  + ^*^3'’’‘2  ' 

(180-4) 

V\'ith  the 

computed 

coefficients  known,  the  optimum  value  of  a 

from  the  relation 

can  be 

J^(a)  = 0 , 

(181) 

which  implies 

that 

"0  = 

(l/lk^)  [-k2+  - (^2-  Sk^k^)]  . 

(182) 

Then,  two  possibilities  arise,  depending  on  the  sign  of 
J^{a)  at  the  point 


(i) 

J (a^)  > 0 , 

fx  0 

(103-1) 

(ii) 

J (n  ) < 0 . 

0 

(183-2) 

In  Case  (i),  the  cubic  interpo] at  ion  process  is  repeated 
between  and  . In  Case  I ii),  it  is  repeated  between  and  1X2  • 
The  process  is  continued  until  Ineq.  (166)  is  satisfied. 
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Limiting  Case.  If  the  solution  of  Eqs.  (180)  is  such 

that 


= 0 , (184) 

the  optimal  stepsize  cannot  be  computed  with  Eq.  (182), 
since  both  the  numerator  and  the  denominator  vanish  simulta- 
neously. This  difficulty  can  be  bypassed  by  observing  that 
the  limiting  case  (184)  means  that  the  cubic  approximation 
(178)  is  being  replaced  by  the  quadratic  approximation  (169). 

As  a consequence,  the  optimal  stepsize  of  the  cubic  approxi- 
mation (182)  must  be  replaced  by  the  optimal  stepsize  of  the 
quadratic  approximation  (174). 

In  practice,  two  cases  arc  possible; 

(i)  ^3  " ^4  ' (185-1) 

(ii)  , (185-2) 

where  i:  ^ is  a small,  preselected  number.  In  Case  (i),  the 
optimal  stepsize  must  be  computed  with  Eq . (182).  In  Case 

(ii),  the  optimal  stepsize  must  bo  computed  with  (174). 

Remark . For  more  details  concerning  the  one-dimensional 
search  for  the  conjugate  gradient  stepsize,  the  reader  is 
referred  to  Ref.  16. 
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9 . Descent  Property  of  a Cycle 

While  the  stepsizes  employed  in  the  conjugate  gradient 
phase  and  the  restoration  phase  are  not  necessarily  small,  a 
descent  property  can  be  proven  for  a complete  conjugate 
gradient-restoration  cycle  under  the  assumption  of  small  step- 
sizes  . 

Let  the  subscript  g denote  the  conjugate  gradient  phase, 
and  let  the  subscript  r denote  the  restoration  phase.  Simple 
manipulations,  omitted  for  the  sake  of  brevity,  show  that  the 
conjugate  gradient  corrections  and  the  restoration  corrections 
have  the  following  orders  of  magnitude: 


Ax^  (t)  = OCotg)  , 

Au  (t) 

g 

= 0(ag), 

All  = 0 (a  ) , 

g g 

(186) 

Ax^  (t)  = , 

Au^(t) 

= 0 (aj-Clg)  - 

2 

Att  = o(a  a ) . 
r I*  g 

(187) 

For  a sufficiently  small,  the  restoration  corrections 
9 

are  negligible  with  respect  to  the  conjugate  gradient  correc- 
tions. Therefore,  providing  the  conjugate  gradient  phase  has 
a descent  property  on  J (a) (this  is  guaranteed  through  the 
selection  of  the  directional  coefficient  y),  the  restoration 
phase  preserves  the  descent  property  of  the  conjugate  gradient 
phase . 

More  specifically,  let  the  subscripts  1,  2,  3 denote 
the  values  of  the  functional  I at  the  beginning  of  the 
conjugate  gradient  phase  , at  the  end  of  the  conju- 
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gate  gradient  phase,  and  at  the  end  of  the  subsequent  restora- 
tion phase.  Vie  note  that  and  I2  are  not  comparable,  since 
the  constraints  are  not  satisfied  to  the  same  accuracy.  On 
the  other  hand,  not  only  and  are  comparable,  but  the 
conjugate  gradient  stepsize  can  be  selected  so  that 

I3  < . (188) 


This  constitutes  the  descent  property  of  a complete  conjugate 
gradient-restoration  cycle. 

If  Ineq.  (188)  is  satisfied,  the  next  conjugate  gradient- 
restoration  cycle  can  be  started.  If  Ineq.  (188)  is  violated, 
one  must  return  to  the  previous  conjugate  gradient  phase  and 
bisect  the  conjugate  gradient  stepsize  until,  after  resto- 
ration, Ineq.  (188)  is  satisfied.  That  the  above  procedure 
leads  to  satisfaction  of  Inec.  (188)  is  guaranteed  by  two 
circumstances:  first,  the  fact  that  the  directional  coefficient 
Y has  been  chosen  consistently  with  Ineq.  (155);  second,  the 
fact  that,  for  small,  the  restoration  corrections  (187)  are 
negligible  by  comparison  with  the  conjugate  gradient  correc- 


tions (186 ) . 
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10 . Summary  of  the  Algorithm 

A sequential  conjugate  gradient-restoration  algorithm 
has  been  developed  in  order  to  solve  optimal  control  problems 
involving  a functional  (12) , subject  to  differential  cons- 
traints (13),  nondif ferontial  constraints  (14),  and  terminal 
constraints (15) . The  algorithm  is  composed  of  a sequence  of 
cycles,  each  cycle  consisting  of  two  phases,  a conjugate 
gradient  phase  and  a restoration  phase. 

Decision  Variables.  The  major  decision  variables  con- 
trolling the  algorithm  are  the  constraint  error  P,  given  by 
Eq.  (24),  and  the  optimality  condition  error  Q,  given  by  Eq. 
(25)  . 

Depending  on  the  value  of  P,  two  cases  are  possible: 


(i) 

P > , 

(189) 

ii) 

P ' C.  ^ 

(190) 

In  Case  (i),  the  algorithm  executes  a restoration  phase.  In 
Case  (ii),  the  algorithm  computes  the  optimality  condition 
error  Q. 

Depending  on  the  value  of  Q,  two  subcases  of  Case  (ii) 


are  possible: 
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(iii)  P < , Q > ^2  , (191 ) 

(iv)  P 


In  Case  (iii),  the  algoritfim  executes  the  conjugate  gradient 
phase.  In  Case  (iv),  the  algorithm  stops:  convergence  has 
been  achieved. 

Iterations . Each  iteration  of  the  conjugate  gradient 
phase  or  the  restoration  phase  is  described  by  the  following 
relations : 


x(t)=x(t)+aA(t),  u(t)=u(t)+-'(D(t),  7i  = TT  + aC,  (193) 


which  tie  the  nominal  functions  and  the  varied  functions. 
Therefore,  each  iteration  includes  two  distinct  operations: 
the  determination  of  the  basic  functions  A(t),  B(t),  C and 

the  determination  of  the  stepsizc  a . 

Restoration  Phase.  The  restoration  phase  includes  one 
or  more  restorative  iterations.  A restorative  iteration  is 
started  whenever  the  constraint  error  P satisfies  Ineq.  (189). 

In  each  restorative  iteration,  the  basic  functions  A(t), 
B(t),  C are  determined  by  solving  the  linear,  two-point 
boundary-value  problem  (43) -(49)  using  the  method  of  particu- 
lar solutions.  This  requires  executing  n+p+1  independent 


sweeps  of  the  system  (43) -(49). 

The  stepsize  a must  be  determined  so  that  the  following 
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1 

inequalities  are  satisfied: 

j P (a)  < P (0)  , T (a)  > 0 . (194) 

For  this  purpose,  a bisection  process,  starting  from  the 
1 reference  stepsize 

I ^0=1,  (195) 

[ is  employed. 

[ 

In  the  course  of  a restorative  iteration,  the  reduction 
of  the  constraint  error  is  guaranteed.  However,  there  is  no 
guarantee  that  the  constraint  error  is  reduced  below  the 
threshold  (190)  characterizing  the  beginning  of  the  next  con- 
jugate gradient  phase.  In  other  words,  after  Ineqs.  (194) 
have  been  satisfied,  two  cases  are  possible: 


(i) 

P (a)  > c , 

(196) 

(ii) 

r (m)  1 . 

(197) 

In  Case  (i),  a further  restorative  iteration  is  initiated 
employing  as  nominal  functions  the  varied  functions  of  the 
previous  restorative  iteration.  In  Case  (ii),  the  restora- 
tion phase  is  terminated,  and  the  next  conjugate  gradient 


phase  is  started. 

Clearly,  each  restoration  phase  includes  a variable 
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number  of  restorative  iterations , depending  on  the  particular 
problem  and  the  nominal  functions  employed.  Generally  speak- 
ing, the  number  of  restorative  iterations  per  restoration 
phase  decreases  in  subsequent  cycles  of  the  sequential  con- 
jugate gradient-restoration  algorithm  and  approaches  zero  as 
the  algorithm  proceeds  toward  convergence. 

Conjugate  Gradient  Phase.  The  conjugate  gradient  phase 
involves  a single  iteration.  This  single  iteration  is  started 
whenever  the  constraint  error  P satisfies  Ineq.  (190). 

In  each  conjugate  gradient  iteration,  the  first  step  is 
to  compute  the  auxiliary  functions  A^{t),  Bj^(t),  C*  corresponding 
to  the  fictitious  value 


Y*=0  (198) 

of  the  directional  coefficient.  These  auxiliary  functions  are 
determined  by  solving  the  linear,  two-point  boundary-value 
problem  (144)-(150)  using  the  method  of  particular  solutions. 
Once  more,  this  requires  executing  n + p+1  independent  sweeps 
of  the  system  (144)-(150). 

With  the  auxiliary  functions  known,  the  basic  functions 
are  determined  with  the  relations 

A(t)=  A*  (t)  + YA(t)  , B(t)=  B*  (t)+  YB(t)  , C = C*+-yC,  (199) 

where  y denotes  the  actual  value  of  the  directional 
coefficient.  This  directional  coefficient  is  sot  at  one  of 
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the  following  levels: 


(i) 

> = 0 , 

(200) 

(ii) 

r = Q/Q  , 

(201) 

where  (200)  holds  for  the  first  conjugate  gradient  iteration 
and  (201)  holds  for  any  subsequent  conjugate  gradient  itera- 
tion . 

Prior  to  accepting  the  directional  coefficient  (200)  or 
(201),  a check  must  be  made.  For  the  choice  (200)  or  (201), 
is  the  slope  of  the  augmented  functional  (0)  negative?  In 
other  words,  does  the  descent  property  of  the  gradient  phase 
hold? 

Concerning  Case  (i),  two  subcases  are  possible: 


(iii) 

1 = 0, 

J^(0)  ^ 0 , 

(202) 

(iv) 

r = 0 , 

J (0)  > 0 . 

CX 

(203) 

In  Case  (iii),  the  directional  coefficient  (200)  is  accepted, 
and  the  algorithm  completes  the  ordinary  gradient  phase.  In 
Case  (iv),  the  descent  property  of  the  ordinary  gradient  phase 
does  not  hold,  and  the  value  of  the  augmented  functional  can- 
not be  reduced,  owing  to  numerical  inaccuracies;  hence,  this 
constitutes  a nonconvergence  condition  for  the  algorithm  as  a 


whole . 
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Concerning  Case  (ii),  two  subcases  are  possible: 


(V) 

Y = Q/Q  , 

J,,(0)  < 0 , 

(204) 

(vi) 

Y = Q/Q  , 

K 

0 

1 

o 

(205) 

In  Case  (v) , the  directional  coefficient  (201)  is  accepted, 

and  the  algorithm  completes  the  conjugate  gradient  phase.  In 

Case  (vi),  the  directional  coefficient  (201)  is  rejected,  and 

is  replaced  by  the  directional  coefficient  (200).  This  means 

that  the  algorithm  is  restarted  with  an  ordinary  gradient 

phase,  characterized  by  > = 0 . 

After  the  directional  coefficient  has  been  selected,  the 

functions  (199)  are  known,  and  the  one-parameter  family  of 

solutions  ( 1 93 ) can  be  formed.  Then,  the  conjugate  gradient 

stepsize  a must  be  determined  through  a one-dimensional  search 

on  the  augmented  functional  J(a)in  such  a way  that  the  fol- 

9 10 

lowing  inequalities  are  satisfied:  ' 

< ^3  (206) 


Note  that  Ineq.  (208)  prevents  the  constraint  error  P from 
becoming  too  largo  during  the  conjugate  gradient  phase. 

^For  the  details  of  the  one-dimensional  search  technique 
leading  to  the  satisfaction  of  (207)-(209),  see  Ref.  16. 
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and 


J(a) 

< J (0)  , 

(207) 

P (a) 

< P*  ' 

(208) 

1 (a) 

> 0 . 

(209) 

For  this  purpose,  the  cubic  interpolation  procedure  of 
Section  8,  is  employed,  with  an  optional  switch  to  quadratic 
interpolation,  if  needed. 

Conjugate  Gradient-Restoration  Cycle.  Generally  speak- 
ing, the  first  cycle  of  the  algorithm  is  a half  cycle,  in 
that  it  includes  a restoration  phase  only.  Every  subsequent 
cycle  is  a complete  cycle,  in  that  it  includes  both  a conju- 
gate gradient  phase  and  a restoration  phase. 

Between  the  endpoints  of  a complete  conjugate  gradient- 
restoration  cycle,  the  following  descent  property  must  be 
satisfied : 

I3  < , (210) 

where  denotes  the  value  of  the  functional  (12)  at  the 
beginning  of  the  cycle  and  denotes  the  value  of  (12)  at 
the  end  of  the  cycle. 

If  Ineq.  (210)  holds,  the  next  conjugate  gradient- 
restoration  cycle  can  be  started.  If  Ineq.  (210)  is  violated. 
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one  must  return  to  the  previous  conjugate  gradient  phase  and 
bisect  the  conjugate  gradient  stepsize  until,  after  restora- 
tion, Ineq.  (210)  is  satisfied. 
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11 . Remarks  and  Safeguards 

In  this  final  section,  we  include  miscellaneous  consi- 
derations, relevant  to  the  computer  implementation  of  the 
sequential  conjugate  gradient-restoration  algorithm.  We 
also  list  some  important  safeguards. 

Nondifferential  Constraint.  For  the  restoration  phase, 
the  linear,  two-point  boundary-value  problem  (43) -(49)  must 
be  solved.  During  the  execution  of  a sweep,  assume  that  time 
station  t has  been  reached  and  that  A(t),  C,  X(t)  are  known  at 
that  time  station.  Then,  Eqs.  (44)  and  (47)  constitute  a 
system  of  m + k equations  in  the  m + k components  of  the  vectors 
B(t),  p(t).  The  system  admits  a unique  solution  providing  the 
following  relation  is  satisfied: 


whore  I denotes  the  m x m identity  matrix  and  0 denotes  the 
k X k null  matrix. 

An  analogous  remak  holds  for  the  conjugate  gradient 
phase:  here,  the  linear,  two-point  boundary-value  problem 
(144) -(150)  must  be  solved.  During  the  execution  of  a sweep, 
assume  that  time  station  t has  been  readied  and  that  Aj^(t),C^  , 
^^(t)are  known  at  that  time  station.  Then,  Lqs.  (145)  and 


(148)  constitute  a system  of  m + k equations  in  the  m+k 
components  of  the  vectors  Bj^(t),  p*(t)  . Once  more,  the  system 
admits  a unique  solution  providing  relation  (211)  is  satisfied. 

The  implication  of  (211)  is  that,  while  the  state  x 
and/or  the  parameter  tt  can  be  absent  from  the  nondifferential 
constraint  (14),  the  control  u can  never  be  absent.  In  fact, 
u must  be  present  in  each  of  the  k scalar  components  of  the 
vector  S.  Therefore,  suitable  transformations  must  be  intro- 
duced in  order  to  convert  problems  where  the  function  S does 
not  involve  the  control  into  problems  where  the  function  S 
involves  the  control.  For  a discussion  of  these  transforma- 
tions, see  Ref.  9. 

Starting  Condition.  The  present  algorithm  can  be  started 
with  nominal  functions  x(t)  u(t),  t satisfying  condition  (15-1) 
and  violating  none,  one,  or  all  of  conditions  (13),  (14), 

(15-2).  If  the  nominal  functions  are  such  that  Ineq.  (189)  is 
satisfied,  the  algorithm  starts  with  a restoration  phase; 
hence,  the  first  cycle  is  a half  cycle,  including  a restoration 
phase  only.  On  the  other  hand,  if  the  nominal  functions  arc 
such  that  Ineq.  (190)  is  satisfied,  the  algorithm  starts 
with  a conjugate  gradient  phase;  hence,  the  first  cycle  is  a 
complete  cycle,  including  both  a conjugate  gradient  phase  and 
a restoration  phase. 
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Bypassing  Condition.  At  the  end  of  the  conjugate 
gradient  phase  of  any  cycle,  the  constraint  error  P must  be 
computed.  If  Ineq.  (189)  is  satisfied,  a restoration  phase 
is  started.  If  Ineq.  (190)  is  satisfied,  the  restoration 
phase  is  bypassed,  and  the  next  cycle  of  the  algorithm  is 
started . 

Convergence  Conditions.  For  the  restoration  phase  taken 
individually,  convergence  is  achieved  whenever  Ineq.  (190)  is 
satisfied.  For  the  sequential  conjugate  gradient-restoration 
algorithm  taken  as  a whole, convergence  is  achieved  whenever 
Ineqs.  (192)  are  satisfied  simultaneously. 

Safeguards . Let  N denote  the  number  of  iterations.  Within 

each  restoration  phase,  let  denote  the  number  of  restorative 

iterations.  Within  each  restorative  iteration,  let  N,  denote 

br 

the  number  of  bisections  of  the  restoration  stepsize  required 
to  satisfy  Ineqs.  (194).  For  the  conjugate  gradient  phase, 
let  denote  the  number  of  bisections  of  the  conjugate  gra- 

dient stepsize  required  to  satisfy  Ineqs.  (207)-(209).  Finally, 
for  a complete  conjugate  gradient-restoration  cycle,  let 
denote  the  number  of  bisections  of  the  conjugate  gradient 
stepsize  required  to  satisfy  the  cycle  descent  property  (210). 

With  this  understand ing , the  following  safeguards  are 
essential  to  the  proper  implementation  of  the  sequential  con- 
jugate gradient-restoration  algorithm: 
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ii) 

N < N*  , 

(212) 

(ii) 

N < N 
r - r* 

(213) 

(iii) 

Nu  '■  N,  , 

br  - br* 

(214) 

(iv) 

bg  - bg* 

(215) 

(V) 

be  - be* 

(216) 

In  the  above  inequalities,  the  right-hand  sides  are  specified 
upper  bounds. 

Restarting  Conditions.  The  directional  coefficient  of 
the  conjugate  gradient  phase  must  be  reset  at  the  level 

V = 0 (217) 

if  any  of  those  circumstances  arise: 

(i) 

(ii) 

(iii) 

(iv) 

(V) 

Satisfaction  of  Ineq.  (218)  indicates  a loss  of  the 


,(o)  :o , 

> = Q/Q  , 

(218) 

'bg  - ^ ' 

Y = 0 or 

II 

o 

o> 

(219) 

A 1 

U 

XI 

> = 0 or 

Y = Q/Q  , 

(220) 

U > N. 

bg  bg* 

Y = Q/Q  , 

(221) 

be  ^bc*  ' 

Y = Q/Q  . 

(222) 

d 
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descent  property  of  the  present  conjugate  gradient  phase. 

Hence,  the  directional  coefficient  must  be  reset  at  the  level 

(217),  characteristic  of  an  ordinary  gradient  phase. 

Satisfaction  of  Ineq.  (219)  indicates  that  the  optimum 

stepsize  of  the  present  conjugate  gradient  phase  [this  is  the 

stepsize  satisfying  Ineq.  (206)]  cannot  be  employed,  owing  to 

violation  of  one  or  more  of  Ineqs.  (207)-(209).  Hence,  this 

stepsize  must  be  bisected  N,  times  so  as  to  arrive  at  satis- 

bg 

faction  of  (207)-(209),  while  violating  (206).  Because  of 
the  ensuing  large  violations  of  the  orthogonality  and  conju- 
gacy  conditions,  the  directional  coefficient  of  the  next 
conjugate  gradient  phase  must  be  reset  at  the  level  (217), 
character isitc  of  an  ordinary  gradient  phase. 

Satisfaction  of  Ineq.  (220)  indicates  that  the  optimum 
stepsize  of  the  present  conjugate  gradient  phase  cannot  be 
employed,  owing  to  violation  of  the  cycle  descent  property 
(210).  Hence,  this  stepsize  must  be  bisected  times,  so  as 

to  arrive  at  satisfaction  of  (210),  while  violating  (206). 

Once  more,  because  of  the  ensuing  large  violations  of  the 
orthogonality  and  conjugacy  conditions,  the  directional 
coefficient  of  the  next  conjugate  gradient  phase  must  be  reset 
at  the  level  (217),  characteristic  of  an  ordinary  gradient 
phase . 


L. 
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Satisfaction  of  Ineq.  (221)  has  a stronger  implication 
than  satisfaction  of  Ineq.  (219).  It  requires  restarting  with 
Y = 0 in  the  present  conjugate  gradient  phase,  rather  than  the 
next  conjugate  gradient  phase. 

Satisfaction  of  Ineq.  (222)  has  a stronger  implication 
than  satisfaction  of  Ineq.  (220).  It  requires  restarting  with 
> = 0 in  the  present  conjugate  gradient  phase,  rather  than  the 
next  conjugate  gradient  phase. 

Nonconvergence  Conditions.  Tne  sequential  conjugate  gradient- 
restoration  algorithm  must  be  programmed  to  stop  whenever  any 
of  several  circumstances  arise: 


(i) 

N - N*  , 

(223) 

(ii) 

N > H 
r r* 

(224) 

(iii) 

br  br* 

(225) 

(iv) 

J^(0)  > 0, 

Y = 0 , 

(226) 

(V) 

bg  bg* 

> = 0 , 

(227) 

(vi) 

be  be* 

Y = 0 , 

(228) 

(vii) 

M > M*  . 

(229) 

Satisfaction  of  Ineq.  (223)  indicates  extreme  slowness  of 
convergence  of  the  algorithm  as  a whole.  Satisfaction  of  Ineq. 


(224)  indicates  extreme  slowness  of  convergence  of  the  resto- 
ration phase.  Satisfaction  of  Ineq.  (225)  indicates  extreme 
smallness  of  the  restorative  displacements.  Satisfaction  of 
Ineq.  (226)  indicates  loss  of  the  descent  property  of  the 
ordinary  gradient  phase,  owing  to  numerical  inaccuracies. 
Satisfaction  of  either  Ineq.  (227)  or  Ineq.  (228)  indicates 
extreme  smallness  of  the  ordinary  gradient  displacements. 
Finally,  satisfaction  of  Ineq.  (229)  indicates  overflow:  the 
modulus  M of  some  of  the  quantities  used  in  the  algorithm  has 
reached  the  upper  limit  M*  allowed  by  the  particular  computer 
employed . 

Rem.rrk.  Several  numerical  examples  illustrating  the 
theory  given  in  this  paper  are  presented  in  Ref.  16.  For  a 
general  discussion  of  the  properties  of  sequential  gradient- 
restoration  algorithms,  the  reader  is  referred  to  Ref.  17. 
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